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APPLICATION OF BOUNDARY INTEGRAL METHOD TO ELASTOPLASTIC 
ANALYSIS OF V-NOTCHED BEAMS 

by Walter Rzasnicki and Alexander Mendelson* 

Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio 

SUMMARY 

The boundary integral equation method was applied in the 
solution of the plane elastoplastic problems. The use of this 
method was illustrated by obtaining stress and strain distri- 
butions for a number of specimens with a single edge notch and 
subjected to pure bending. The boundary integral equation 
method reduced the non-homogeneous biharmonic equation to two 
coupled Fredholm-type integral equations. These integral equations 
were replaced by a system of simultaneous algebraic equations and 
solved numerically in conjunction with the method of successive 
elastic solutions. 

INTRODUCTION 

Knowledge of the stress distribution in the neighborhood 
of a singularity, such as the tip of a V-notch in a bar loaded 
in tension or bending, is of fundamental importance in evaluating 
the resistance to fracture of structural materials. Elastic 
solutions to various geometries have been obtained by a number 
of different methods. Among the more effective ones, are the 

complex variable method (ref. 1), collocation method (ref. 2), 

_ 
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and finite element method (ref. 3) . However, the first two of 
these methods are not general enough nor readily adaptable to 
three-dimensional or elastoplastic problems. And the finite 
element method requires solutions of large sets of equations and 
fails to give sufficiently fine resolution in the vicinity of 
notch tips . 

The recently developed boundary integral methods (ref. 4) 
offer an attractive alternative to other methods of analysis. 

These methods have a number of advantages as listed in references 4 
and 5, the most important of which is that nodal points are needed 
only on the boundary instead of throughout the region as required 
by finite element methods. To date these methods have been used 
primarily for obtaining elastic solutions to various problems. 

Their extension to elastoplastic problems has been proposed in 
references 4 and 6. However no elastoplastic solution has here- 
tofore been obtained. 

The present paper describes the application of a boundary in- 
tegral method to the solution of the elastic and elastoplastic 
problems of a V-notched beam in pure bending. Solutions of such 
problems for the elastoplastic case by finite elements has not 
been too successful in obtaining fine enough resolution and 
sufficiently accurate results in the vicinity of the notch tip 
(ref. 7). The present method overcomes these difficulties, 
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since the stress and strain can be computed at any arbitrary point 
in the body from a knowledge of boundary values only. 

METHOD OF SOLUTION 

The problem of determining the state of stress and strain 
in a plane elastoplastic problem can be reduced to solving the 
following non-homogeneous biharmonic equation for the Airy stress 
function, as shown in reference 8: 

V 4 <P = g (x,y) (1) 


g(x,y) = - 


l - y 
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for the plane strain case, and 


g(x,y) = - E 


— r- (E P + Ae P ) + — j (E P + Ae P ) 
3y 2 X X 3x 2 y y 
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3x3y 
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for the plane stress case, where e P , e P , and e P y represent the 
accumulation of plastic strain increments from the beginning 
of the loading history up to, but not including the current 



4 
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increment of the load, and Ae , Ae„, and Ae r are the increments 

x Y xy 

of plastic strain due to the current increment of load. For the 
elastic case, g(x,y) is of course equal to zero. 

The stress function cp must satisfy appropriate boundary 
conditions. For the problem under consideration (fig. 1) <P (x,y) 
and its outward normal derivative 3cp/3n must satisfy the fol- 
lowing boundary conditions (ref. 9) : 


along boundary OA and OA' <p(x,y) = 0; — = 0 

3n 




dtp 

along boundary AB and A'B' cp(x,y) =0; — = 0 

3n 


along boundary BC and B'C' 


/( 4 ) 
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along boundary CD and C'D cp(x,y) = ; — = 0 

3n 


6 



5 


To solve equation (1) by means of the boundary integral 
method use is made of Green's second theorem to reduce this 
equation to coupled integral equations, as shown in references 4 
and 9. The result is 
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where 


4 = V cp 


p = r In r 


and 

r(x,y;£,n) is the distance between any two points P(x,y) and 
qU,n) or P (x,y) and QU,n), as shown in figure 2. 


( 8 ) 


Equation (5) would, for a known function g(x,y), give us 

directly a solution to the biharmonic equation (1) provided the 

2 2 

functions cp(x,y) , 9cp(x,y)/9n, V cp(x,y) and 9[V <P(x,y)]/9n 
were known on the boundary C. 

However, only the stress function cp and its outward normal 

derivative 9<P/9n are specified (eq. (4)). The values of 
2 2 

V cp = <t> and 9 (V cp)/3n = 9<J>/9n on the boundary must be com- 
patible with the given values of <P and 9<P/3n. To assure 
this compatibility, we have to solve the system of coupled 
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integral equations (6) and (8) , which contain the unknown 
functions $ and 3$/3n. 

Once the values of $ and 3$/3n on the boundary C of region R 
are known we can proceed with the calculation of the stress field in 
the region R utilizing equation (5) and the equations which define cp » 
namely 


a 

x 


2 2 
3 cp/3y , 


2 2 

a-3 cp/3x o Y . r 
y xy 


2 

3 cp/3x3y 


1 9 ) 


The calculation of the function g(x,y) which is obtained 
iterc- ively, is discussed in detail in reference 9. 

NUMERICAL PROCEDURES 
Solution of the Integral Equations 
Since it is generally impossible to solve the system of 
coupled integral equations analytically, a numerical method is 
utilized in which the integral equations (6) and (8) are re- 
placed by a system of simultaneous algebraic equations. 

For simplicity of notation the normal derivatives are 
denoted by prime superscripts. The boundary is divided into n 
intervals, not necessarily equal, numbered consecutively in the 
direction of increasing s. The center of each interval is de- 
signated as a node. The values of $ and 4> ' are assumed con- 
stants on each interval and equal to the values calculated at the 


node . 
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In a similar manner the interior of region R is covered by a 
grid, containing m cells. The cells do not have to have equal areas. 
Their nodal points are located at the centroids. The value of 
g(5,n) is assumed constant over each cell and equal to the value 
calculated at the centroid. 

The arrangements of boundary and interior region sub-divisxons 
is shown in figures 3 and 4. Using the above assumptions, 
equations (6) and (8) can be replaced by a system of 2n simulta- 


neous algebraic equations with 2n unknowns, that is, 4> and <t>^ 
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where i = 1,2,3. ..n, r.. is the distance from i fc ^ node to the 

xk 

centroid of the k th cell, is the area of the k th cell, 
and 




a. . = 
ID 


(In r . . ) ' ds 
ID 


In r. .ds 
ID 


p . .ds 
ID 


d. . = - / p . .ds 

ID . / ID 


e. . = / (V p. .) 'ds 

iD ./ ID 


f . . = - / V p. .ds 

ID / ID 


where integration is taken over the j interval, and r^ is the 
distance from i th node to any point in the j th interval. The 
normal derivatives in equations (11) are taken on the j th interval. 

For curved boundaries the coefficients given by equations (11) 
can be evaluated, if necessary, by Simpson's rule for i i- j. For 
i = j, because of the singular nature of the integrand, the inte- 
grals for the coefficients must be evaluated by a limiting process. 


For boundary intervals, such as for the problem treated herein, 
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which can be represented by straight lines a closed form solution can 
be obtained for these coefficients (ref. 9) . 

Boundary equations (10) can be written in matrix form as 

[B] {x} = {r} (12) 

where [B] is 2n x 2n matrix and (x} and are 2n x 1 column 
matrices . 

Matrix [B] is dependent only on geometry, that is, number of 
nodes and their distribution on the boundary. Since the matrix ( r ) 
contains the nonlinear function g(£,n), which depends on the stress 
field and therefore on matrix {x}, and iterative process will be 
used to obtain the solution. 

To calculate stresses, at any nodal point in the region R, 

from the stress function cp we need not perform any numerical 
differentiation. Equation (5) can be differentiated under the 
integral sign and once $ and $ ' are known on the boundary 
the stresses can be obtained by the same type of numerical inte- 
gration as in equations (10) . Applying equations (9) to equation (5) 
yields for the case of a rectangular grid the following stress 


equations : 
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where now i = L,2,3...m and refers to the centroid of the i 

cell, 6^ and 6^ represent, respectively, x-directional and 

y-directional dimension of che cell. The coefficients A. . , 

B.., C . . , D.., E.., F.., G.., H.,, I.., K.. are obtained hy 
il il i] il il il il il il 

appropriate differentiation under the integral sign of the 
coefficients given by equations (11) and are listed in 
reference 9. 

Boundary Interval and Interior Grid Size 
Use was made of the geometric and loading symmetry about 
the x axis for the problem created herein to reduce the number of 
unknowns by approximately a factor of two. Since the vicinity 
of the notch tip is of greatest interest, a fine nodal spacing 
along the notch was chosen. However to maintain a reasonable 
number of nodal points, and at the same time to obtain fine re- 
solution at the tip of the notch, the boundary along the notch 
was divided into a number of intervals progressively increasing 
in length as one moved away from the tip. The rate of change 
in th interval's length along this boundary was optimized 
by the method presented in reference 5. For the cases con- 
sidered optimum ratios of the lengths of two consecutive boundary 
intervals were found to be in the range of 1.08 to 1.10. The 
resulting smallest dimensionless boundary interval length varied 
from 0.0001 to 0.0002. The nodal arrangement shown in figure 5 
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was used for all cases considered, resulting in a set of 140 
equations containing 140 unknowns. 

The choice of the size of the grid, which has to cover the 
region where plastic flow is expected to occur, is discussed 
in detail in reference 9. The smallest cells located in the 
vicinity of the tip of the notch have dimensions 6^/w = 0.004; 

/w = 0.008 for plane strain cases, and 5 A*' = 0.004, 

y c x 

- /w = 0.016 for Diane stress cases A tyoical interi.r ^ria 

y 

is shown in figure 5. 

The solution to this boundary value problem is obtained 
by an iterative process, known as the method of successive 
elastic solutions described in detail in reference 8. This 
method, applied to the present problem, is presented in ref- 
erence 9. Once the successive approximation procedure has con- 
verged, the stresses and strains are known everywhere in the beam. 

At,SULTS AND DISCUSSION 

A number of beam problems were solved for both plane stress 
and plane strain cases. These included notch depth to beam depth 
ratios of 0.3 and 0.5, notch angles of 3° and 10°, strain hard- 
ening parameter values of 0.05 and 0.1. In addition calculations 
were performed using the actual stress-strain curve of a 5083-0 
aluminum alloy in order to compare the calculated results for 
the crack opening displacements with available experimental 
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results. The detailed results for all cases are given in 
reference 9. Here, only a few typical results are presented. 

Calculations were first made to check the method for the 
elastic problem, with the function g set equal to zero. For 
the elastic case a larger number of notch depths was used than 
for the elastoplastic case. A comparison with the boundary 
collocation results of reference 2 is shown in figures 6 to 8. 

Very good agreement was obtained. 

The stress intensity factor Kj is here defined by 

K_ = iim r n o (r,0) (14) 

L r-0 y 

Figures 9 to 12 show typical plastic zone growth with load 
for a plane strain case and for a plane stress case. 

Variation of the dimensionless stress intensity factor with 
load is shown in figure 13 for the case of a specimen with notch 
depth of a/w = 0.5 anc c = 10°, under plane strain condition 
and two values of strain hardening parameter m. It is to be noted 
that n appearing in equation (14) is no longer constant as in the 
elastic case, but is a function of the applied load. The stress 
intensity factor shows no significant increase over the linear 
elastic value up to an applied load of q=o,/o rt =0.40. 

Above this load increases progressively for both m's, at the 

faster rate for lower strain hardening parameter. 
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In order to verify in part, the numerically computed 
results, the notch opening displacements for a specimen with 
a 10° edge notch, notch depth a/w =0.5 and a stress- strain 
curve given by figure 14 were compared with experimental re- 
sults obtained by Bubsey, R. T. and Jones, M. (ref. 10). The 
specimen used in this experiment, made of A1 5083-0 with length 
to width ratio of 4:1, crack length a/w = 0.5 was subjected 
to three point bending. The experimental results as shown in 
figure 15 are in good agreement with the numerical results ob- 
tained herein. 

The product of y-directional stress and total strain 
was calculated for various cases. The order of singularity of 
that product was determined by plotting lnCa^e^) versus In r 
and by making a least squares fit of a straight line through the 
plotted points. It was found to be very close to unity. 

Finally, Rice's J integral was evaluated for several 
cases by using relations given in reference 5. In these calcu- 
lations, straight line paths were chosen through the plastic zone 
near the tip of the notch. The integral was evaluated using 
values of stresses, strains and displacements at cells centroids 
for a number of paths. The path- independence of J was not con- 
clusive, since the results varied up to 15 percent from the averaged 
value. It is possible that the results obtained herein do not 
indicate that the path independent property is lost but rather that 
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the field values of the displacements are not calculated with 
sufficient accuracy. 

The average values of the dimensionless J integral as a 
function of load are plotted in figure 16 for a case of a speci- 
menn with a 10° edge notch, a/w = 0.5, m = 0.05 and plane strain 
condition. At the start of plastic flow J increases rapidly 
with load. This is followed by almost linear variation with 
additional “load. 

The relations between the J integral and stress inten- 
sity factor Kj developed for linear elasticity are obviously not 

applicable for the elastoplastic problem. By plotting the 
2 

dimensionless J/K^. ratios as a funtion of load q, the relation 

between the J integral and the dimensionless stress intensity 

factor K J is obtained. A typical plot is shown in figure 17. 

2 

In all cases, the ratio J/K ^ remains almost equal to elastic 
value of 0.89 for plane strain or 1.0 for plane stress and in- 
creases sharply at the load corresponding to the appearance 
of the plastic zone at the boundary opposite the notch. Once 
the transition occurs the ratio increases approximately propor- 
tionally to the load increment. 
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CONCLUSIONS 

The boundary integral equation method proved to be ca- 
pable of giving very detailed results such as stress and strain 
distributions around the tip of the notch and, related to them, 
the shapes of plastic zones. This was accomplished using a rela- 
tively ,mall number unknowns. 

The obtained rerults also provide the information of 
the effect of strain hardening and on the differences that occur 
between plane stress and plane strain solutions. The order of 
singularity for the strain energy density was found to be unity, 
which is consistent with results previously obtained by other 
investigators. 

The stress intensity factor was calculated for several 
cases. The path independence property of Rice's J integral 
was qualitatively confirmed and the relation between J and the 
stress intensity factor was graphically extended into the elasto- 
plastic range. 

The presence of a singularity at the tip of the notch 
makes accurate answers very difficult to obtain. Nevertheless 
good agreement was obtained between the calculated results and 
experimentally measured notch opening displacement as shown in 
figure 15. Some improvement in the solution techniques and 
further investigation of the influence of the boundary nodal 
spacing and interior grid size on the resulting stress and 
strain fields, and therefore, on the notch opening displacements 
and J integrals, may be desirable. 
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Figure 16. - Dimensionless plane strain J Integral tor specimen with a lOf 3 edge 
notch subjected to pure bending; a/w ■ 0. 5, a • 10P, m • 0. 05, p ■ 0, 31 



fioure 17. - Vacation ol the ratio ol dimensionless J Integral to the square ol 
dimensionless stiess intensity factor with load lor a specimen with a 10P edge 
notch subjected to pure bending; plane strain, a/w • 0. 5, a • loP, m * 0, 35, p * 0. 31 
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